TP 2 - Extensions du modèle linéaire
library(bestglm)
library(ggfortify)
library(ggplot2)
library(gridExtra)
library(VGAM)
library(ROCR)
library(tidyverse)
library(leaps)
library(MASS)
library(glmnet)
library(plotly)
library(corrplot)Sélection de variables
Les données
On reprend le jeu de données “ozone.txt” du TP1.
Les 13 variables observées sont :
- maxO3 : Maximum de concentration d’ozone observé sur la journée en \(\mu\)gr/m3
- T9, T12, T15 : Température observée à 9, 12 et 15h
- Ne9, Ne12, Ne15 : Nébulosité observée à 9, 12 et 15h
- Vx9, Vx12, Vx15 : Composante E-O du vent à 9, 12 et 15h
- maxO3v : Teneur maximum en ozone observée la veille
- vent : orientation du vent à 12h
- pluie : occurrence ou non de précipitations
Les données sont disponibles ici.
Pour les lire, vous pouvez utiliser la commande suivante :
url <- url("https://raw.githubusercontent.com/cmaugis/Stat-Avancees---PFB-Toulouse/main/Data/Ozone.txt")
ozone <- read.table(url)
ozone$vent <- as.factor(ozone$vent)
ozone$pluie <- as.factor(ozone$pluie)Vous pouvez obtenir un résumé du jeu de données ozone avec la commande summary().
maxO3 T9 T12 T15
Min. : 42.00 Min. :11.30 Min. :14.00 Min. :14.90
1st Qu.: 70.75 1st Qu.:16.20 1st Qu.:18.60 1st Qu.:19.27
Median : 81.50 Median :17.80 Median :20.55 Median :22.05
Mean : 90.30 Mean :18.36 Mean :21.53 Mean :22.63
3rd Qu.:106.00 3rd Qu.:19.93 3rd Qu.:23.55 3rd Qu.:25.40
Ne9 Ne12 Ne15 Vx9
Min. :0.000 Min. :0.000 Min. :0.00 Min. :-7.8785
1st Qu.:3.000 1st Qu.:4.000 1st Qu.:3.00 1st Qu.:-3.2765
Median :6.000 Median :5.000 Median :5.00 Median :-0.8660
Mean :4.929 Mean :5.018 Mean :4.83 Mean :-1.2143
3rd Qu.:7.000 3rd Qu.:7.000 3rd Qu.:7.00 3rd Qu.: 0.6946
Vx12 Vx15 maxO3v vent pluie
Min. :-7.878 Min. :-9.000 Min. : 42.00 Est :10 Pluie:43
1st Qu.:-3.565 1st Qu.:-3.939 1st Qu.: 71.00 Nord :31 Sec :69
Median :-1.879 Median :-1.550 Median : 82.50 Ouest:50
Mean :-1.611 Mean :-1.691 Mean : 90.57 Sud :21
3rd Qu.: 0.000 3rd Qu.: 0.000 3rd Qu.:106.00
[ getOption("max.print") est atteint -- ligne 1 omise ]On réajuste ici un modèle de régression linéaire avec toutes les variables explicatives (cf TP1) pour la suite.
Call:
lm(formula = maxO3 ~ ., data = oz_quanti)
Residuals:
Min 1Q Median 3Q Max
-53.566 -8.727 -0.403 7.599 39.458
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.24442 13.47190 0.909 0.3656
T9 -0.01901 1.12515 -0.017 0.9866
T12 2.22115 1.43294 1.550 0.1243
T15 0.55853 1.14464 0.488 0.6266
Ne9 -2.18909 0.93824 -2.333 0.0216 *
Ne12 -0.42102 1.36766 -0.308 0.7588
Ne15 0.18373 1.00279 0.183 0.8550
Vx9 0.94791 0.91228 1.039 0.3013
Vx12 0.03120 1.05523 0.030 0.9765
Vx15 0.41859 0.91568 0.457 0.6486
maxO3v 0.35198 0.06289 5.597 1.88e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.36 on 101 degrees of freedom
Multiple R-squared: 0.7638, Adjusted R-squared: 0.7405
F-statistic: 32.67 on 10 and 101 DF, p-value: < 2.2e-16Question 1
Enoncé
Dans le summary(reg.mul), un test est fait sur chaque coefficient. Il revient à tester que la variable n’apporte pas d’information supplémentaire sachant que toutes les autres variables sont dans le modèle.
Ici, nous souhaitons aller plus loin dans la possible simplification du modèle. Nous allons utiliser des procédures de choix de modèles pour sélectionner les variables significatives. On va ici comparer la sélection de variable obtenue par différents critères: BIC, AIC, \(R^2\) ajusté, \(C_p\) de Mallows. Pour cela, vous pouvez utiliser la fonction regsubsets() de la librairie leaps et la fonction stepAIC() de la librairie MASS. Commentez les résultats obtenus avec les différents critères, vous pourrez vous aider des commandes suivantes :
Correction
choix <- regsubsets(maxO3 ~ ., data = oz_quanti, nbest = 1, nvmax = 11, method = "backward")
plot(choix, scale = "Cp")
Start: AIC=607.26
maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 +
maxO3v
Df Sum of Sq RSS AIC
- T9 1 0.1 20827 605.26
- Vx12 1 0.2 20827 605.26
- Ne15 1 6.9 20834 605.30
- Ne12 1 19.5 20847 605.36
- Vx15 1 43.1 20870 605.49
- T15 1 49.1 20876 605.52
- Vx9 1 222.6 21050 606.45
<none> 20827 607.26
- T12 1 495.5 21323 607.89
- Ne9 1 1122.6 21950 611.14
- maxO3v 1 6459.6 27287 635.51
Step: AIC=605.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Vx12 1 0.1 20827 603.26
- Ne15 1 6.9 20834 603.30
- Ne12 1 22.1 20849 603.38
- Vx15 1 43.5 20871 603.49
- T15 1 49.4 20877 603.52
- Vx9 1 264.0 21091 604.67
<none> 20827 605.26
- T12 1 641.4 21469 606.66
- Ne9 1 1183.6 22011 609.45
- maxO3v 1 7112.2 27940 636.16
Step: AIC=603.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Ne15 1 6.9 20834 601.30
- Ne12 1 23.1 20850 601.38
- T15 1 50.0 20877 601.53
- Vx15 1 79.6 20907 601.69
- Vx9 1 320.1 21148 602.97
<none> 20827 603.26
- T12 1 647.3 21475 604.69
- Ne9 1 1185.1 22013 607.46
- maxO3v 1 7112.6 27940 634.16
Step: AIC=601.3
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Ne12 1 16.2 20851 599.38
- T15 1 45.2 20880 599.54
- Vx15 1 75.1 20910 599.70
- Vx9 1 325.2 21160 601.03
<none> 20834 601.30
- T12 1 971.8 21806 604.40
- Ne9 1 1215.8 22050 605.65
- maxO3v 1 7106.3 27941 632.17
Step: AIC=599.38
maxO3 ~ T12 + T15 + Ne9 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- T15 1 45.6 20896 597.63
- Vx15 1 77.6 20928 597.80
- Vx9 1 337.6 21188 599.18
<none> 20851 599.38
- T12 1 1034.0 21885 602.80
- Ne9 1 2155.5 23006 608.40
- maxO3v 1 7130.2 27981 630.33
Step: AIC=597.63
maxO3 ~ T12 + Ne9 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Vx15 1 74.1 20970 596.02
- Vx9 1 366.5 21263 597.58
<none> 20896 597.63
- Ne9 1 2241.1 23137 607.04
- T12 1 6712.6 27609 626.83
- maxO3v 1 7381.5 28278 629.51
Step: AIC=596.02
maxO3 ~ T12 + Ne9 + Vx9 + maxO3v
Df Sum of Sq RSS AIC
<none> 20970 596.02
- Vx9 1 903.4 21874 598.75
- Ne9 1 2714.8 23685 607.66
- T12 1 6650.4 27621 624.88
- maxO3v 1 7363.5 28334 627.73
Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = oz_quanti)
Coefficients:
(Intercept) T12 Ne9 Vx9 maxO3v
12.6313 2.7641 -2.5154 1.2929 0.3548 Avec le \(C_p\) de Mallows, le critère BIC et le R2 ajusté, on retient 4 variables explicatives : T12, Ne9, Vx9 et maxO3v.
Même concluson en utilisant la fonction stepAIC.
On retrouve une cohérence avec les corrélations entre variables
Question 2
Enoncé
A la question précédente, nous retenons 4 variables explicatives T12, Ne9, Vx9 et maxO3v.
A l’aide des commandes suivantes, testez le sous-modèle avec les 4 variables retenues contre le modèle complet. Commentez.
Correction
Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = oz_quanti)
Residuals:
Min 1Q Median 3Q Max
-52.396 -8.377 -1.086 7.951 40.933
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.63131 11.00088 1.148 0.253443
T12 2.76409 0.47450 5.825 6.07e-08 ***
Ne9 -2.51540 0.67585 -3.722 0.000317 ***
Vx9 1.29286 0.60218 2.147 0.034055 *
maxO3v 0.35483 0.05789 6.130 1.50e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14 on 107 degrees of freedom
Multiple R-squared: 0.7622, Adjusted R-squared: 0.7533
F-statistic: 85.75 on 4 and 107 DF, p-value: < 2.2e-16Analysis of Variance Table
Model 1: maxO3 ~ T12 + Ne9 + Vx9 + maxO3v
Model 2: maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 +
maxO3v
Res.Df RSS Df Sum of Sq F Pr(>F)
1 107 20970
2 101 20827 6 143.01 0.1156 0.9944La p-valeur est \(p=0.994 \geq 0.05\) donc on ne rejette pas l’hypothèse nulle au risque \(5\%\) donc on retient le modèle simplifié avec 4 variables explicatives.
Régressions régularisées
On va s’intéresser aux régressions ridge, Lasso et Elastic-net dans cette partie. On commence par centrer et réduire les données avant de mettre en place et comparer ces méthodes de régression régularisées.
tildeY <- scale(oz_quanti[, 1], center = TRUE, scale = TRUE)
tildeX <- scale(oz_quanti[, -1], center = TRUE, scale = TRUE)Régression Ridge
Question 1
Enoncé
A l’aide de la fonction glmnet(), ajustez une régression ridge en faisant varier \(\tau\) sur une grille. On stockera le résultat dans la variable fitridge. Explorez le contenu de fitridge.
Correction
tau_seq <- seq(0, 1, by = 0.001)
fitridge <- glmnet(tildeX, tildeY, alpha = 0, lambda = tau_seq, family = "gaussian",
intercept = FALSE)
summary(fitridge) Length Class Mode
a0 1001 -none- numeric
beta 10010 dgCMatrix S4
df 1001 -none- numeric
dim 2 -none- numeric
lambda 1001 -none- numeric
dev.ratio 1001 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 7 -none- call
nobs 1 -none- numericfitridge$lambda contient les valeurs du paramètre \(\tau\) dans l’ordre décroissant. fitridge$beta contient la valeur des paramètres \(\theta_j\) pour chaque valeur de \(\tau\). Comme on a centré les données, on n’a pas d’intercept donc fitridge$a0 est un vecteur de zéro.
Pour l’interprétation des autres sorties, vous pouvez vous reporter à l’aide de la fonction glmnet.
Question 2
Enoncé
Tracez les chemins de régularisation de chaque variable et commentez.
Question 3
Enoncé
A l’aide de la fonction cv.glmnet() mettez en place une validation croisée pour sélectionner le “meilleur” \(\tau\) par l’erreur quadratique moyenne (MSE).
df2 <- data.frame(tau = ridge_cv$lambda, MSE = ridge_cv$cvm, cvup = ridge_cv$cvup,
cvlo = ridge_cv$cvlo)
ggplot(df2) + geom_line(aes(x = tau, y = MSE)) + geom_vline(xintercept = ridge_cv$lambda.min,
col = "red", linetype = "dotted") + geom_line(aes(x = tau, y = cvup), col = "blue",
linetype = "dotted") + geom_line(aes(x = tau, y = cvlo), col = "blue", linetype = "dotted") +
xlim(c(0, 1))La valeur de \(\tau\) sélectionnée vaut ….
Correction
ridge_cv <- cv.glmnet(tildeX, tildeY, alpha = 0, lambda = tau_seq, nfolds = 10, type.measure = "mse",
intercept = FALSE)
best_tau <- ridge_cv$lambda.min
best_tau[1] 0.239df2 <- data.frame(tau = ridge_cv$lambda, MSE = ridge_cv$cvm, cvup = ridge_cv$cvup,
cvlo = ridge_cv$cvlo)
ggplot(df2) + geom_line(aes(x = tau, y = MSE)) + geom_vline(xintercept = ridge_cv$lambda.min,
col = "red", linetype = "dotted") + geom_line(aes(x = tau, y = cvup), col = "blue",
linetype = "dotted") + geom_line(aes(x = tau, y = cvlo), col = "blue", linetype = "dotted") +
xlim(c(0, 1))
La valeur de \(\tau\) sélectionnée vaut 0.239.
g1 <- g1 + geom_vline(xintercept = best_tau, linetype = "dotted", color = "red") +
xlim(c(0, best_tau + 0.1))
g1
Régression Lasso
Question 1
Enoncé
A l’aide de la fonction glmnet(), ajustez une régression Lasso en faisant varier \(\tau\) sur une grille. On stockera le résultat dans la variable fitlasso. Explorez le contenu de fitlasso.
Correction
tau_seq <- seq(0, 1, 0.001)
fitlasso <- glmnet(tildeX, tildeY, alpha = 1, lambda = tau_seq, family = "gaussian",
intercept = FALSE)
summary(fitlasso) Length Class Mode
a0 1001 -none- numeric
beta 10010 dgCMatrix S4
df 1001 -none- numeric
dim 2 -none- numeric
lambda 1001 -none- numeric
dev.ratio 1001 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 7 -none- call
nobs 1 -none- numericVoici quelques exemples de valeurs estimées pour les paramètres pour plusieurs valeurs de \(\tau\)
Question 2
Enoncé
Tracez le chemin de régularisation de chacune des variables et commentez.
Question 3
Enoncé
A l’aide de la fonction cv.glmnet() mettez en place une validation croisée pour sélectionner le “meilleur” \(\tau\) par MSE.
lasso_cv <- cv.glmnet(...) # A COMPLETER
best_tau <- lasso_cv$lambda.min
best_tau1se <- lasso_cv$lambda.1seLa valeur de \(\tau\) sélectionnée vaut ….
Correction
La valeur de \(\tau\) sélectionnée vaut 0.026 ( 0.164 pour un écart d’un s.e du min)
g2 <- g2 + geom_vline(xintercept = best_tau, linetype = "dotted", color = "red") +
geom_vline(xintercept = best_tau1se, linetype = "dotted", color = "green")
g2
Question 4
Correction
11 x 1 sparse Matrix of class "dgCMatrix"
s1
(Intercept) .
T9 .
T12 0.300068858
T15 0.079165203
Ne9 -0.144351908
Ne12 -0.002279795
Ne15 .
Vx9 0.039391227
Vx12 .
Vx15 .
maxO3v 0.251543269Régression Elastic Net
Question 1
Enoncé
Reprenez les questions précédentes pour ajuster une régression Elastic Net
fitEN <- glmnet(...)
df <- data.frame(tau = rep(fitEN$lambda, ncol(tildeX)), theta = as.vector(t(fitEN$beta)),
variable = rep(c(colnames(tildeX)), each = length(fitEN$lambda)))
g3 <- ggplot(df, aes(x = tau, y = theta, col = variable)) + geom_line() + theme(legend.position = "bottom")EN_cv <- cv.glmnet(tildeX, tildeY, alpha = 0.5, lambda = tau_seq, nfolds = 10, type.measure = "mse",
intercept = FALSE)
g3 <- g3 + geom_vline(xintercept = EN_cv$lambda.min, linetype = "dotted", color = "red") +
geom_vline(xintercept = EN_cv$lambda.1se, linetype = "dotted", color = "green")
ggplotly(g3)Correction
tau_seq <- seq(0, 1, 0.001)
fitEN <- glmnet(tildeX, tildeY, alpha = 0.5, lambda = tau_seq, family = "gaussian",
intercept = FALSE)
df <- data.frame(tau = rep(fitEN$lambda, ncol(tildeX)), theta = as.vector(t(fitEN$beta)),
variable = rep(c(colnames(tildeX)), each = length(fitEN$lambda)))
g3 <- ggplot(df, aes(x = tau, y = theta, col = variable)) + geom_line() + theme(legend.position = "bottom")EN_cv <- cv.glmnet(tildeX, tildeY, alpha = 0.5, lambda = tau_seq, nfolds = 10, type.measure = "mse",
intercept = FALSE)
g3 <- g3 + geom_vline(xintercept = EN_cv$lambda.min, linetype = "dotted", color = "red") +
geom_vline(xintercept = EN_cv$lambda.1se, linetype = "dotted", color = "green")
ggplotly(g3)Question 2
Enoncé
Comparez les coefficients obtenus avec la régression linéaire, la régression ridge, la régression Lasso et la régression ElasticNet
Correction
Modèle linéaire généralisé
Régression logistique
Dans cette section, nous allons nous intéresser au jeu de données SAheart disponible dans la librairie bestglm. On souhaite expliquer la présence/absence d’une maladie cardiovasculaire (chd) en fonction de 9 variables et on dispose pour cela d’un échantillon de \(n=462\) individus.
Question 1
A l’aide des commandes suivantes, chargez les données et faites quelques statistiques descriptives pour vous familiariser avec les données.
La description des données est disponible dans l’aide de R (tapez ?SAheart dans la console).
Le jeu de données est composé des variables suivantes :
- sbp: pression artérielle systolique
- tabac: tabac cumulatif (kg)
- ldl: cholestérol des lipoprotéines de basse densité
- adiposité:
- famhist: antécédents familiaux de maladie cardiaque (Présent, Absent)
- type: comportement de type-A
- obésité:
- alcool: consommation actuelle d’alcool
- âge: âge au début
- chd: réponse, maladie coronarienne
'data.frame': 462 obs. of 10 variables:
$ sbp : int 160 144 118 170 134 132 142 114 114 132 ...
$ tobacco : num 12 0.01 0.08 7.5 13.6 6.2 4.05 4.08 0 0 ...
$ ldl : num 5.73 4.41 3.48 6.41 3.5 6.47 3.38 4.59 3.83 5.8 ...
$ adiposity: num 23.1 28.6 32.3 38 27.8 ...
$ famhist : Factor w/ 2 levels "Absent","Present": 2 1 2 2 2 2 1 2 2 2 ...
$ typea : int 49 55 52 51 60 62 59 62 49 69 ...
$ obesity : num 25.3 28.9 29.1 32 26 ...
$ alcohol : num 97.2 2.06 3.81 24.26 57.34 ...
$ age : int 52 63 46 58 49 45 38 58 29 53 ...
$ chd : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 1 2 1 2 ... sbp tobacco ldl adiposity
Min. :101.0 Min. : 0.0000 Min. : 0.980 Min. : 6.74
1st Qu.:124.0 1st Qu.: 0.0525 1st Qu.: 3.283 1st Qu.:19.77
Median :134.0 Median : 2.0000 Median : 4.340 Median :26.11
Mean :138.3 Mean : 3.6356 Mean : 4.740 Mean :25.41
3rd Qu.:148.0 3rd Qu.: 5.5000 3rd Qu.: 5.790 3rd Qu.:31.23
Max. :218.0 Max. :31.2000 Max. :15.330 Max. :42.49
famhist typea obesity alcohol age
Absent :270 Min. :13.0 Min. :14.70 Min. : 0.00 Min. :15.00
Present:192 1st Qu.:47.0 1st Qu.:22.98 1st Qu.: 0.51 1st Qu.:31.00
Median :53.0 Median :25.80 Median : 7.51 Median :45.00
Mean :53.1 Mean :26.04 Mean : 17.04 Mean :42.82
3rd Qu.:60.0 3rd Qu.:28.50 3rd Qu.: 23.89 3rd Qu.:55.00
Max. :78.0 Max. :46.58 Max. :147.19 Max. :64.00
chd
0:302
1:160
SAheart_quanti <- SAheart[, -c(5, 10)]
A <- data.frame(y = as.vector(as.matrix(t(SAheart_quanti))), x = rep(colnames(SAheart_quanti),
nrow(SAheart_quanti)))
g1 <- ggplot(A, aes(x = x, y = y)) + geom_boxplot()
g2 <- ggplot(SAheart, aes(x = chd)) + geom_bar()
g3 = ggplot(SAheart, aes(x = famhist)) + geom_bar()
print(g1)
Question 2
Enoncé
Ajustez un modèle de régression logistique multiple additif à l’aide de la fonction glm().
Correction
Modèle de régression logistique :
\[chd_i\sim\mathcal{B}(\pi_i)\] avec
\[g(\pi_i)=\ln\left(\frac{\pi_i}{1-\pi_i}\right) = \theta_0 + \sum_{j=1}^{9} \theta_j x_i^{(j)} \textrm{ dont } x_i^{(famhist)} = \mathbb{1}_{famhist_i = present}\]
Call:
glm(formula = chd ~ ., family = binomial(link = "logit"), data = SAheart)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7781 -0.8213 -0.4387 0.8889 2.5435
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.1507209 1.3082600 -4.701 2.58e-06 ***
sbp 0.0065040 0.0057304 1.135 0.256374
tobacco 0.0793764 0.0266028 2.984 0.002847 **
ldl 0.1739239 0.0596617 2.915 0.003555 **
adiposity 0.0185866 0.0292894 0.635 0.525700
famhistPresent 0.9253704 0.2278940 4.061 4.90e-05 ***
typea 0.0395950 0.0123202 3.214 0.001310 **
obesity -0.0629099 0.0442477 -1.422 0.155095
alcohol 0.0001217 0.0044832 0.027 0.978350
age 0.0452253 0.0121298 3.728 0.000193 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 596.11 on 461 degrees of freedom
Residual deviance: 472.14 on 452 degrees of freedom
AIC: 492.14
Number of Fisher Scoring iterations: 5Question 3
Correction
[1] 0.2079628La valeur du pseudo-\(R^2\) n’est pas très élevé.
Question 4
Correction
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.1507208650 1.308260018 -4.70145138 2.583188e-06
sbp 0.0065040171 0.005730398 1.13500273 2.563742e-01
tobacco 0.0793764457 0.026602843 2.98375801 2.847319e-03
ldl 0.1739238981 0.059661738 2.91516648 3.554989e-03
adiposity 0.0185865682 0.029289409 0.63458325 5.257003e-01
famhistPresent 0.9253704194 0.227894010 4.06052980 4.896149e-05
typea 0.0395950250 0.012320227 3.21382267 1.309805e-03
obesity -0.0629098693 0.044247743 -1.42176449 1.550946e-01
alcohol 0.0001216624 0.004483218 0.02713729 9.783502e-01
age 0.0452253496 0.012129752 3.72846442 1.926501e-04On a en dernière colonne les pvaleurs associées à chaque \(Z\)-test pour tester la nullité de chacun des coefficients. On peut remarquer que les variables sbp, adiposity, obesity et alcohol ont une pvaleur assez élevée. On va poursuivre à l’étude en cherchant à simplifier le modèle par sélection de variables.
Question 5
Enoncé
Comparez les différentes procédures de sélection de variables mises en place par les commandes suivantes:
Correction
On commence par une procédure de sélection de variables avec la fonction bestglm() en utilisant le critère BIC :
Call: glm(formula = y ~ ., family = family, data = Xi, weights = weights)
Coefficients:
(Intercept) tobacco ldl famhistPresent typea
-6.44644 0.08038 0.16199 0.90818 0.03712
age
0.05046
Degrees of Freedom: 461 Total (i.e. Null); 456 Residual
Null Deviance: 596.1
Residual Deviance: 475.7 AIC: 487.7Fitting algorithm: BIC-glm
Best Model:
df deviance
Null Model 456 475.6856
Full Model 461 596.1084
likelihood-ratio test - GLM
data: H0: Null Model vs. H1: Best Fit BIC-glm
X = 120.42, df = 5, p-value < 2.2e-16Le modèle retenu est \[ logit(\pi_i) = \theta_0 + \theta_2\ tobacco_i + \theta_3\ ldl_i + \theta_5\ \mathbb{1}_{famhist_i = present} + \theta_6\ typea_i + \theta_9\ age_i \]
On fait de même avec le critère AIC :
Call: glm(formula = y ~ ., family = family, data = Xi, weights = weights)
Coefficients:
(Intercept) tobacco ldl famhistPresent typea
-6.44644 0.08038 0.16199 0.90818 0.03712
age
0.05046
Degrees of Freedom: 461 Total (i.e. Null); 456 Residual
Null Deviance: 596.1
Residual Deviance: 475.7 AIC: 487.7On sélectionne le même modèle.
On peut aussi faire la procédure de sélection de variables à l’aide de la fonction step(). Par défaut, on utilise le critère AIC. On peut voir les étapes de la procédure de sélection backward. On retrouve le même sous-modèle retenu.
Start: AIC=492.14
chd ~ sbp + tobacco + ldl + adiposity + famhist + typea + obesity +
alcohol + age
Df Deviance AIC
- alcohol 1 472.14 490.14
- adiposity 1 472.55 490.55
- sbp 1 473.44 491.44
<none> 472.14 492.14
- obesity 1 474.23 492.23
- ldl 1 481.07 499.07
- tobacco 1 481.67 499.67
- typea 1 483.05 501.05
- age 1 486.53 504.53
- famhist 1 488.89 506.89
Step: AIC=490.14
chd ~ sbp + tobacco + ldl + adiposity + famhist + typea + obesity +
age
Df Deviance AIC
- adiposity 1 472.55 488.55
- sbp 1 473.47 489.47
<none> 472.14 490.14
- obesity 1 474.24 490.24
- ldl 1 481.15 497.15
- tobacco 1 482.06 498.06
- typea 1 483.06 499.06
- age 1 486.64 502.64
- famhist 1 488.99 504.99
Step: AIC=488.55
chd ~ sbp + tobacco + ldl + famhist + typea + obesity + age
Df Deviance AIC
- sbp 1 473.98 487.98
<none> 472.55 488.55
- obesity 1 474.65 488.65
- tobacco 1 482.54 496.54
- ldl 1 482.95 496.95
- typea 1 483.19 497.19
- famhist 1 489.38 503.38
- age 1 495.48 509.48
Step: AIC=487.98
chd ~ tobacco + ldl + famhist + typea + obesity + age
Df Deviance AIC
- obesity 1 475.69 487.69
<none> 473.98 487.98
- tobacco 1 484.18 496.18
- typea 1 484.30 496.30
- ldl 1 484.53 496.53
- famhist 1 490.58 502.58
- age 1 502.11 514.11
Step: AIC=487.69
chd ~ tobacco + ldl + famhist + typea + age
Df Deviance AIC
<none> 475.69 487.69
- ldl 1 484.71 494.71
- typea 1 485.44 495.44
- tobacco 1 486.03 496.03
- famhist 1 492.09 502.09
- age 1 502.38 512.38Pour utiliser le critère BIC avec la fonction step(), il faut utiliser la syntaxe suivante :
Question 6
Enoncé
Confirmez par un test de sous-modèle le sous-modèle retenu modelbest dans la question précédente.
Par quels critères peut-on aussi comparer les deux modèles modelbest et modellogit ?
Correction
On commence par ajuster le modèle retenu par les méthodes de sélection de variables
modelbest <- glm(chd ~ tobacco + ldl + famhist + typea + age, family = binomial,
data = SAheart)
summary(modelbest)
Call:
glm(formula = chd ~ tobacco + ldl + famhist + typea + age, family = binomial,
data = SAheart)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9165 -0.8054 -0.4430 0.9329 2.6139
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.44644 0.92087 -7.000 2.55e-12 ***
tobacco 0.08038 0.02588 3.106 0.00190 **
ldl 0.16199 0.05497 2.947 0.00321 **
famhistPresent 0.90818 0.22576 4.023 5.75e-05 ***
typea 0.03712 0.01217 3.051 0.00228 **
age 0.05046 0.01021 4.944 7.65e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 596.11 on 461 degrees of freedom
Residual deviance: 475.69 on 456 degrees of freedom
AIC: 487.69
Number of Fisher Scoring iterations: 5On peut par exemple remarquer que le test de nullité de chacun des coefficients séparement a une p-valeur \(\ll 0.05\).
On fait ensuite un test de sous-modèle entre modelbest et modellogit
Analysis of Deviance Table
Model 1: chd ~ tobacco + ldl + famhist + typea + age
Model 2: chd ~ sbp + tobacco + ldl + adiposity + famhist + typea + obesity +
alcohol + age
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 456 475.69
2 452 472.14 4 3.5455 0.471La pvaleur valant 0.47, on ne rejette pas le sous-modèle au risque \(5\%\). Pour le sous-modèle on a \(k_0=6\) paramètres à estimer donc \(n-k_0=462-6=456\)
Pour le modèle complet on a \(k_1=10\) paramètres à estimer donc \(n-k_0=462-10=452\)
\(\mathcal{D}(M_0) - \mathcal{D}(M_1) \stackrel{\mathcal{L}}{\longrightarrow} \chi^2(k_1-k_0 =4)\)
Critère pseudo-\(R^2\) :
[1] 0.202015[1] 0.2079628Critère AIC :
[1] 487.6856[1] 492.14Question 7
Enoncé
Comment s’interprètent les différents coefficients de modelbest ?
(Intercept) tobacco ldl famhistPresent typea
0.001586152 1.083693731 1.175850406 2.479793436 1.037812583
age
1.051755195 Correction
Pour l’interprétation des paramètres, il faut considérer des odds-ratio, attention entre variables qualitatives et quantitatives.
Régression loglinéaire
Dans cette partie, on souhaite étudier la diversité des fourmis sur le site expérimental des Nourragues en Guyane Française. Les données, disponibles dans le fichier Fourmis.txt, sont issues du protocole expérimental suivant : des morceaux de 1\(m^2\) de litière ont été récoltés. Ils ont été pesés (car le poids de litière est vu comme un indicateur de l’épaisseur de la litière) et le nombre d’espèces différentes présentes dans l’échantillon a été dénombré. Enfin les conditions de recueil (humides ou sèches) ont été notées pour tester leur influence sur la présence des fourmis. L’objectif est de mettre en évidence les variables qui influencent potentiellement le nombre d’espèces de fourmis présentes dans le milieu.
Question 1
Chargez le jeu de données Fourmis et executez les commandes suivantes pour faire quelques statistiques descritpives.
url <- url("https://raw.githubusercontent.com/cmaugis/Stat-Avancees---PFB-Toulouse/main/Data/Fourmis.txt")
Fourmis <- read.table(url, header = TRUE, sep = ",")
Fourmis$Site <- as.factor(Fourmis$Site)
Fourmis$Conditions <- as.factor(Fourmis$Conditions)
summary(Fourmis) Effectifs Weight Site Conditions
Min. : 2.00 Min. : 6.200 FLWT:50 Dry:83
1st Qu.: 7.00 1st Qu.: 8.525 FTWT:50 Wet:87
Median :10.00 Median : 9.950 GPWT:50
Mean :11.19 Mean :10.168 INWT:20
3rd Qu.:15.00 3rd Qu.:11.600
Max. :28.00 Max. :17.800 # A tibble: 8 x 3
# Groups: Site [4]
Site Conditions `n()`
<fct> <fct> <int>
1 FLWT Dry 24
2 FLWT Wet 26
3 FTWT Dry 24
4 FTWT Wet 26
5 GPWT Dry 22
6 GPWT Wet 28
7 INWT Dry 13
8 INWT Wet 7
Question 2
Enoncé
Dans cette question, on cherche à expliquer le nombre de fourmis présentes dans un échantillon de sol en prenant en compte les conditions de recueil, le site et l’interaction entre ces deux facteurs.
- Ecrivez un modèle linéaire généralisé modpois adapté et programmez-le à l’aide de la fonction
glm().
Etudiez si on peut simplifier le modèle modpois. Vous pouvez vous aider de la fonction
anova(....,test="Chisq").A l’aide des commandes suivantes, on obtient le nombre moyen d’espèces de fourmis attendu sur un échantillon de terre aux différents sites pour les deux conditions. Commentez.
Correction
On considère un modèle loglinéaire modpois pour expliquer le nombre de fourmis présentes dans un échantillon de sol en prenant en compte les conditions de recueil, le site et l’interaction entre ces deux facteurs.
\[ Eff_i\sim\mathcal{P}(\lambda_i) \textrm{ et } \ln(\lambda_i) = \theta_0 + \sum_{j\in S}(\theta_j + \gamma_j \mathbb{1}_{Conditions_i=WET})\mathbb{1}_{Site_i=j} + \alpha \mathbb{1}_{Conditions_i=WET} \] avec \(S=\{FTWT,GPWT,INTW\}\).
Call:
glm(formula = Effectifs ~ Site * Conditions, family = poisson,
data = Fourmis)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.7219 -1.2101 -0.2415 0.9066 2.9152
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.82583 0.04969 56.869 < 2e-16 ***
SiteFTWT -0.62861 0.08425 -7.461 8.60e-14 ***
SiteGPWT -0.69113 0.08857 -7.803 6.06e-15 ***
SiteINWT -0.42794 0.09727 -4.399 1.09e-05 ***
ConditionsWet -0.13851 0.07132 -1.942 0.0521 .
SiteFTWT:ConditionsWet -0.17515 0.12476 -1.404 0.1603
SiteGPWT:ConditionsWet 0.28473 0.11880 2.397 0.0165 *
SiteINWT:ConditionsWet 0.63099 0.14148 4.460 8.20e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 530.36 on 169 degrees of freedom
Residual deviance: 329.68 on 162 degrees of freedom
AIC: 1045.4
Number of Fisher Scoring iterations: 4On étudie si on peut simplifier le modèle modpois. On commence par tester si on peut enlever l’interaction.
anova(glm(Effectifs ~ Site + Conditions, data = Fourmis, family = poisson), modpois,
test = "Chisq")Analysis of Deviance Table
Model 1: Effectifs ~ Site + Conditions
Model 2: Effectifs ~ Site * Conditions
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 165 361.05
2 162 329.68 3 31.367 7.114e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1On rejette l’hypothèse nulle, on conserve l’effet d’interaction. On ne peut donc pas simplifier le modèle.
Nombre moyen d’espèces de fourmis attendu sur un échantillon de terre aux différents sites pour les deux conditions :
newx <- data.frame(Site = rep(levels(Fourmis$Site), 2), Conditions = rep(levels(Fourmis$Conditions),
each = 4))
newx <- data.frame(newx, lambdahat = predict(modpois, newdata = newx, type = "response"))
newx Site Conditions lambdahat
1 FLWT Dry 16.875000
2 FTWT Dry 9.000000
3 GPWT Dry 8.454545
4 INWT Dry 11.000000
5 FLWT Wet 14.692308
6 FTWT Wet 6.576923
7 GPWT Wet 9.785714
8 INWT Wet 18.000000Question 3
Enoncé
On souhaite maintenant étudier le nombre d’espèces de fourmis présentes dans une unité de volume, en fonction des caractéristiques du site. Nous allons ici intégrer l’information du poids (Weight) qui est vu comme un indicateur de l’épaisseur de la litière.
On remarque que la relation entre les paramètres des lois de Poisson du nombre d’espèces présentes en un site et du nombre d’espèces présentes dans une unité de volume est
\[ \mathbb{E}[\textrm{nb esp.}]=\mathbb{E}[\textrm{nb esp. par unité de volume}]\times Weight \]
On en déduit donc un modèle pour le nombre d’espèces par unité de volume :
\[Eff_i\sim\mathcal{P}(\lambda_i)\] avec
\[\ln(\lambda_i) = \ln(Weight_i) + \theta_0 + \sum_{j\in S}(\theta_j + \gamma_j \mathbb{1}_{Conditions_i=WET})\mathbb{1}_{Site_i=j} + \alpha \mathbb{1}_{Conditions_i=WET} \]
La variable Weight est donc considérée comme un offset.
Ajustez ce modèle à l’aide de l’option
"offset"de la fonctionglm().Quel est le nombre moyen d’espèces de fourmis attendu sur un échantillon de terre qui pèse 10kg aux différents sites pour les deux conditions ?
Correction
Modèle pour le nombre d’espèces par unité de volume : \[Eff_i\sim\mathcal{P}(\lambda_i)\] avec
\[\ln(\lambda_i) = \ln(Weight_i) + \theta_0 + \sum_{j\in S}(\theta_j + \gamma_j \mathbb{1}_{Conditions_i=WET})\mathbb{1}_{Site_i=j} + \alpha \mathbb{1}_{Conditions_i=WET} \]
glmInt <- glm(Effectifs ~ Site * Conditions, offset = log(Weight), family = "poisson",
data = Fourmis)
summary(glmInt)
Call:
glm(formula = Effectifs ~ Site * Conditions, family = "poisson",
data = Fourmis, offset = log(Weight))
Deviance Residuals:
Min 1Q Median 3Q Max
-3.2267 -0.9268 -0.1522 0.8252 3.6257
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.49531 0.04969 9.968 < 2e-16 ***
SiteFTWT -0.55200 0.08425 -6.552 5.69e-11 ***
SiteGPWT -0.74225 0.08857 -8.380 < 2e-16 ***
SiteINWT -0.38293 0.09727 -3.937 8.26e-05 ***
ConditionsWet -0.02802 0.07132 -0.393 0.69438
SiteFTWT:ConditionsWet -0.37426 0.12476 -3.000 0.00270 **
SiteGPWT:ConditionsWet 0.16746 0.11880 1.410 0.15866
SiteINWT:ConditionsWet 0.47387 0.14148 3.349 0.00081 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 534.13 on 169 degrees of freedom
Residual deviance: 301.41 on 162 degrees of freedom
AIC: 1017.2
Number of Fisher Scoring iterations: 4Le pseudo-R\(^2\) de ce modèle :
[1] 0.4357083Le nombre moyen d’espèces de fourmis attendu sur un échantillon de terre qui pèse 10kg aux différents sites pour les deux conditions ?
newx <- data.frame(Site = rep(levels(Fourmis$Site), 2), Conditions = rep(levels(Fourmis$Conditions),
each = 4), Weight = rep(10, 8))
newx <- data.frame(newx, lambdahat = predict(glmInt, newdata = newx, type = "response"))
newx Site Conditions Weight lambdahat
1 FLWT Dry 10 16.410049
2 FTWT Dry 10 9.448819
3 GPWT Dry 10 7.811844
4 INWT Dry 10 11.189358
5 FLWT Wet 10 15.956558
6 FTWT Wet 10 6.319290
7 GPWT Wet 10 8.980662
8 INWT Wet 10 17.475728Question 4
Enoncé
Est-ce-qu’une modélisation avec un modèle prenant en compte la sur-dispersion est plus appropriée ?
Correction
modquasipoisson <- glm(Effectifs ~ Site * Conditions, offset = log(Weight), family = quasipoisson(link = "log"),
data = Fourmis)
summary(modquasipoisson)
Call:
glm(formula = Effectifs ~ Site * Conditions, family = quasipoisson(link = "log"),
data = Fourmis, offset = log(Weight))
Deviance Residuals:
Min 1Q Median 3Q Max
-3.2267 -0.9268 -0.1522 0.8252 3.6257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.49531 0.06739 7.350 9.23e-12 ***
SiteFTWT -0.55200 0.11426 -4.831 3.12e-06 ***
SiteGPWT -0.74225 0.12012 -6.179 5.02e-09 ***
SiteINWT -0.38293 0.13191 -2.903 0.00421 **
ConditionsWet -0.02802 0.09672 -0.290 0.77239
SiteFTWT:ConditionsWet -0.37426 0.16918 -2.212 0.02836 *
SiteGPWT:ConditionsWet 0.16746 0.16110 1.039 0.30015
SiteINWT:ConditionsWet 0.47387 0.19186 2.470 0.01455 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 1.839056)
Null deviance: 534.13 on 169 degrees of freedom
Residual deviance: 301.41 on 162 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4Pseudo-R\(^2\) pour ce modèle :
[1] 0.4357083